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Abstract. We present a statistical model which is able to capture some interesting features exhibited in 
the Brazilian test of rock samples. The model is based on elements which break irreversibly when the force 
experienced by the elements exceed their own load capacity. If an element breaks the load capacity of the 
neighboring elements are decreased by a certain amount, assuming weakening effect around the defected 
zone. From the model we numerically investigate the stress-strain behavior, the strength of the system, 
how it scales with the system size and also it’s fluctuation for both uniform and Weibull distribution of 
breaking thresholds in the system. To check the validity of our statistical model we perform few Brazilian 
tests on Sandstone and Chalk samples. The stress-strain curve from model results agree qualitatively well 
with the lab-test data. Also, the damage profile right at the point when the stress-strain curve reaches it’s 
maximum is seen to mimic the crack patterns observed in our Brazilian test experiments. 


1 Introduction 


The tensile strength of the rock is one of the important pa¬ 
rameters that influences the rock crushing and rock blast¬ 
ing results. To measure the tensile strength of a solid body 
(rock, concrete, ceramics etc.) the Brazilian test [lj[2j[3j|4] 
is performed. This test is also named as the Diametral 
Compression test and the Indirect Tensile test. The test 
was introduced by Carniero [5] in Brazil, and Akazawa [6] 
in Japan in 1943. 

In the Brazilian test a circular disk is diametrically 
compressed. Due to this compression a tensile stress de¬ 
velops perpendicular to the direction of the applied force 
and the amplitude of the induced stress is proportional to 
the applied force. When the induced tensile stress crosses 
certain limit, fracture develops mainly in the middle zone 
of the sample. The indirect tensile strength of a cylindrical 
sample with diameter D and thickness t is given by pQ , 


2 F 

at ~ nDt 


( 1 ) 


where F is the force at the failure point of the material. 
This equation was obtained analytically based on the as¬ 
sumption that the rock is isotropic and homogeneous [7, 
[8]. However, the nature is more complex. Therefore we 
must take into account the anisotropy and heterogeneity 
effects for better estimation of rock-strength. 
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Fig. 1 . The loading apparatus used in our experiment. 


Besides the tensile strength, the crack initiation and 
crack propagation upon the Brazilian test is often exam¬ 
ined mm- A lot of studies have been done on the crack 
initiation position in the Brazilian test both experimen¬ 
tally and numerically 031 . In this paper, we present our 
experimental study of Brazilian test on two samples and 
then present a very simple model in the framework of sta¬ 
tistical physics, which is able to capture some features 
exhibited in the Brazilian test. 

In Sect. [2] we describe the experimental procedure of 
our Brazilian test followed by some experimental results 
in Sect. [3j We elaborately describe our model in Sect. 0] 
and present numerical results obtained from this model in 
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Fig. 2. The load-time behavior of Castlegate sandstone (CG) 
and Mons chalk (MC) samples. 

Sect. [5] Here we investigate stress-strain behavior, dam¬ 
age profile and critical strength of the system. Finally we 
discuss the application of this work and make some con¬ 
cluding remarks in Sect. [6] 


2 Experimental Procedure 

The Brazilian (splitting) indirect tensile strength test [8| 
has been performed on two rock types - Castlegate sand¬ 
stone (CG) and Mons chalk (MC). Both rocks have a 
fairly homogeneous and isotropic (with no distinct lay¬ 
ering) structure. The specimens were cored perpendicular 
to bedding and cut into the discs of approximately 22 mm 
thickness and 52 mm diameter. Each disc was wrapped 
around the circumference with a masking tape. This soft¬ 
ened the contact between the rock and the steel curved 
jaws (Fig. |T]) within which the disc is clamped during the 
test. This curved jaws is then placed in a loading frame. 
We then continuously increase the externally applied load, 
at a very slow but constant speed of crosshead, until the 
failure of the sample occurs. In our experiment the ex¬ 
ternal force is applied with a MTS load frame at the rate 
0.003 mm/Sec., until the failure of the material. The reac¬ 
tion force given by the sample is recorded at every time. 
Fracture development, indicated by the sudden drop in 
loading, is captured using a digital camera. 



Fig. 3. Fracture patterns after failure. Top one is for Castle¬ 
gate sandstone (CG) and the lower one is for Mons chalk (MC) 
samples. 

in the Table [lj In our experiment we increase the exter¬ 
nally applied load on the material and measure the force 
with time. In Fig.[2]the load-time data has been plotted for 
both samples. The force/load values are measured in units 
of Newton (N). We define the tensile strength of material 
as the maximum load that a material can sustain with¬ 
out failure. For CG sample the first drop in the measured 


Table 1. Details about the samples that we examined. The 
thickness/diameter ratio is 0.42 for both of the samples. 


3 Experimental Results 

We experimentally examined the strength of two cylindri¬ 
cal samples CG and MC -masses and dimensions are given 


Sample name 

Weight 

Diameter 

Thickness 

Castlegate sandstone 

91.61gm. 

52.26mm. 

21.78mm. 

Mons chalk 

71.73gm. 

51.75mm. 

21.73mm. 
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Fig. 4. (a) Square lattice of size L = 3. Each point on the 
lattice is an element, follows Hooke’s law before the breaking 
point. A force Fo is applied at the middle column from the 
top of the lattice, (b) Illustrates the effect of the applied force. 
Here it is assumed that the effect decreases linearly from the 
point of loading. The cut off length L c and the To determines 
the slope of the decreasing force profile, slope m=-Fo/L c . 

load at around 195 Sec., signifies the failure. When a frac¬ 
ture opens the near-by sand grains occupy the fractured 
region, also the contact area to the curved jaws broadens 
and thus the reaction force increases further. This is not 
observed in case of chalk sample due to the brittleness of 
the chalk. From the load-time plot it is obvious that the 
tensile strength of CG is greater than that of MC. In Fig. 
[3] the damage zone has been shown for both of the samples 
after the failure. 


4 Model 

Our model consists of Hookean elements, placed on each 
lattice point of a regular square lattice, with their own 
breaking thresholds drawn randomly from a given proba¬ 
bility distribution. Once the load on an element exceeds 
its breaking threshold, it breaks immediately and the load 
capacity of its neighboring intact elements are decreased 
indicating the weakening effect around the defected zone. 
One may think of the redistribution of load carried by the 
broken element among the intact elements as in the case 
of fiber bundle model mug. In our model the load expe¬ 
rienced by the intact elements are kept constant through¬ 
out the breaking process, elements break only due to the 




Fig. 5. (a) For the uniformly distributed thresholds between 
0.1 to 1 and L c = L/12 , the stress-strain behavior of a sys¬ 
tem with L — 64 has been shown for three different values of 
strength decrement factor a = 0.50 (red), a = 0.70 (black) and 
a = 0.90 (blue), (b) The same data in (a) are plotted against 
the scaled variable so -0 ' 92 and we obtain a reasonably good 
data collapse upto the failure point (the peak position). 

decrement of their breaking thresholds. To model the Brazil¬ 
ian test, the external load is applied at the middle point 
of the system and it is assumed that the effect of the ap¬ 
plied force diminishes linearly as we move-on horizontally 
in either sides from the point of loading. From this model 
we numerically investigate the breaking properties of the 
system. 

We consider a square lattice of size L as shown in the 
Fig. Ufa). Each point on the lattice is an element. Each 
individual element positioned at (x, y) has been assigned 
a breaking threshold b( x ^ i.e., it can sustain a maximum 
6( Xj2/ ) amount of load through it, beyond which it breaks ir¬ 
reversibly. Before breakdown each element follows Hooke’s 
law. For simplicity we assume the stiffness factor of all the 
elements to be unity. The breaking thresholds {b( XiV )} are 
drawn from a probability distribution p(b), the cumula¬ 
tive distribution of which is given by P(b) = f 0 p(y)dy. 
The breaking threshold distribution is the only source of 
disorder present in this model. 

In our model we apply a force at the middle point of 
the lattice from the top of it, as shown in the Fig.0Ja). We 
assume that the effect of this applied force extends within 
a certain region from the point of application of the force. 
We consider the effect of this force decreases linearly and 
this force profile is symmetric about the point of loading. 
The force profile at the positive part of the lattice has the 
form: 


F(x) = mx + Fq 


( 2 ) 
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Fig. 6. Crack profile for different values of the strength decre¬ 
ment factor, a , right at the point where the stress-strain curve 
reached it’s maximum. The black color represents the intact 
elements and the white patches at the middle of it are the 
broken elements within the system of size L — 64. Here the 
thresholds are distributed uniformly between 0.1 to 1 and in 
simulation we used L c = L/ 12. 


applied force at the point x = 0 and m defines how far the 
effect of this applied force should extend. In our model we 
prefix a cut-off length L c beyond which the applied force 
has no effect i.e., F(x ) = 0 for x > L c . This actually 
defines the value of m(=-Fo/L c ). The whole schematic 
picture is depicted in the Fig. 4(b). We also assume that 
the load through all the elements of a particular column 
is same i.e., F(x,y) = F{x). 

An element breaks irreversibly and it’s neighborhood 
gets affected, when the force through the element exceeds 
it’s breaking threshold. When an element breaks, the break¬ 
ing strength of it’s surviving nearest neighbors is decreased 
by a factor <a(0 < a < 1), i.e., their new thresholds are 
assigned a value equal to a times of their own previous 
thresholds. This is actually the sense of weakening around 
the defected zone. The force profile remains unaltered in 
this model. 

When a force is applied to the system, the breakdown 
of the elements occur from the middle part (close to the 
x = 0) of the system since the force field is much stronger 
in this region. So, the local failures will grow from the 
middle zone of the system. This local failures decrease 
the strength of the elements around the defected zone. A 
local failure may occurs due to the weakness of the ele¬ 
ments as well as due to the decrement of the strength of 
the elements after an failure event and due to this a spa¬ 
tial correlation develops in the failure events. Since the 
force field decreases linearly in the horizontal direction, 
it is more probable that failure will occur in the vertical 
direction. This makes cracks to propagate in the vertical 
up or vertical down direction rather than horizontal di¬ 
rection from the defected zone, which correlates with the 
experimental observations. 

The width of disorder, the value of the cut-off length 
L c and the strength decrement factor a play a crucial role 
in the crack propagation. For a fixed window of disorder, 
smaller the value of a and L c , lead successive failure to 
be more and more localized and directed in the vertical 
direction. 


5 Numerical Results 


We investigate the stress-strain behavior of the system 
which is the most crucial part of the Brazilian test. We 
observe the stress-strain behavior, the system size depen¬ 
dence of the maximum value of the stress and it’s fluc¬ 
tuation for both uniform and Weibull distribution of the 
breaking thresholds. 

To obtain the stress-strain curve we increase the exter¬ 
nal load quasi-statically. Due to the loading all elements 
get stretched. At an elongation, 5 per intact element, the 
total force at a column is e times the total number of 
intact elements at the column. We define stress (a) as- 


(7 = 



(3) 


and since it is symmetric about the point of loading so we where, Ni n is the number of intact elements at the central 
have F(-x) = F(x). In Eqn. (2) Fq is the amplitude of the column. 
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Fig. 7. (a) Plot of cr c (L) with L. The system size varies from 
L — 32 to L — 1024 with an increment of factor 2. (b) Plot of 
cr c (L) with L~ oa is best fitted by a straight line, the extrapo¬ 
lation of which gives the value of a c = 0.14. 


Numerically, the stress-strain curve is obtained in the 
following way. Initially we start with a fully intact sys¬ 
tem and then we apply a load such that the weakest el¬ 
ement in the system breaks. The weakest surviving ele¬ 
ment is determined by calculating the maximum of the 
ratio between force field to the breaking threshold i.e., 
max{F(x,y)/b( Xj y)}, which is equal to max{F(x)/b( x ^}. 
Once it is removed, the breaking strength of it’s nearest 
surviving neighbor is reduced by a factor a. Due to this, it 
may so happen that the load through some of the elements 
exceed their breaking threshold and we break them simul¬ 
taneously again reducing the strength of their neighboring 
elements, that may cause further breakdown and so on. 
This process stops when all the remaining intact elements 
have their breaking strength greater than their respective 
experienced stress values. We then increase the externally 
applied load in such a way that the weakest element among 
the remaining intact system reaches it’s breaking thresh¬ 
old and then the entire procedure is repeated. This pro¬ 
cedure is continued until we are at a point well above the 
point where the maximum of the stress-strain curve takes 
place. 

In Fig. (5ja) the stress-strain curve is shown for a par¬ 
ticular configuration for different values of a. In the next 
figure by suitably scale the abscissa and ordinate we could 
obtain a data collapse (Fig. [5jb)) upto the peak point. In 
experiments, after peak position, different rock samples 
behave differently -depending on their failure structure 
(clear fracture, complex fracture, crashed etc). Therefore, 
in our numerical study, we have tried the scaling up to 
the peak position. The defected zone right at the point 



Fig. 8. Log-log plot of W(L) with L for system size upto 
L a= 512. The fitted line has a slope 0.46 and thus we conclude 
W(L) ~ L~ 0A6 . 

where the stress-strain curve reaches it’s maximum has 
been shown by the white color in Fig. [6j Throughout the 
simulation we have used L c = L/ 12. 

Once we generate the stress-strain curve, we then cal¬ 
culate the critical strength of the system. We define critical 
strength crf(L) for a particular sample s of size L with a 
particular set of breaking thresholds {&( Xj2/ )} as the max¬ 
imum value of the stress in the stress-strain curve. For 
each sample s we numerically calculate cr?(L) and when 
the whole procedure is repeated over a large number of un¬ 
correlated samples, it eventually gives the average value 
of the critical strength a c (L) = ( <J S C (L )) for the system 
of size L. We then repeat the entire set of calculation for 
different values of L to know how cr c (L) approaches to it’s 
asymptotic value. 

We assume that the average value cr c (L) for a system 
of size L converges to it’s asymptotic value a c = a c (oo) as 
L oo according to the following form, 

a c (L) =a c + AL~ X I V (4) 

For uniformly distributed breaking threshold of the ele¬ 
ments in the range [0.1-1.0], the variation of cr c (L) against 
L has been shown in Fig. [7fa). When the same cr c (L) val¬ 
ues are plotted in Fig. Ulb) as a function of L~ x / V with 
1/v = 0.40, we observe that it can be fitted by a straight 
line. The least-squares fitted straight line cuts the vertical 
axis and has the form cr c (L ) = 0.14 + 0.28L -1 / 1 '. From 
here our conclusion is that a c = 0.14 and v = 5/2 for the 
uniformly distributed breaking thresholds. 

We also measure the fluctuation of the critical strength, 
defined as- 

W(L) = vV(L)2) - <<r c (L)>2 (5) 

In Fig. [8] we plot W(L) for different values of L in log-log 
scale for the uniformly distributed breaking threshold and 
we observe that it decreases with increase in L. The best 
fitted line has a slope 0.46. 

We then repeat our entire set of studies for elements 
with Weibull breaking threshold distribution of the form: 

P(b) = 1 - e~ b ' 


(6) 








6 


Sumanta Kundu et al.: A simple Discrete-Element-Model of Brazilian Test 



i- 1 -1- 1 -1- 1 -1- 1 -1- 1 -1- 1 - 

0 0.05 0.1 0.15 0.2 0.25 0.3 



Fig. 9. (a) Plot of cr c (L) with L~ 0A . The best fitted straight 
line when extrapolated, it cuts the vertical axis at 0.24 and 
thus the value of cr c = 0.24. Here the breaking threshold of 
the elements are drawn from the Weibull distribution of the 
form of Eqn. (6). (b) The fluctuation W(L) varies with L as- 
W(L) ~ L -0 ' 51 . 

We have estimated the values of cr c (L) numerically for 
six different system sizes from 32 to 1024 increased by a 
factor of 2 at every step. Plotting them against L~ 0A and 
extrapolating as L —>• oo we have obtained a c = 0.24. In 
Fig. [U(a) the plot of a c (L) with L~ 0A has been shown. Our 
conclusion is that for the Weibull distribution of breaking 
thresholds of the form of Eqn. (6), the value of a c = 0.24 
and v = 5/2. 

We also calculate W(L) for different values of L upto 
L = 512. In Fig. [9](b) we plot the W{L) as a function of 
L and this is fitted with a straight line, resulting W(L) ~ 
L-°- 51 . 


6 Discussions and Conclusions 

To model the Brazilian test we have introduced a statisti¬ 
cal model and studied numerically the strength of the ma¬ 
terial. The stress-strain behavior of the system is found to 
have a maxima. We have investigated the critical strength 
of the system and also the system size dependence of 
the critical strength for both uniformly and Weibull dis¬ 
tributed breaking threshold of the elements in the system. 
We observe that the damage profile right at the critical 
point (Fig. [6j) resembles the nature of the damage profile 
observed in the Brazilian test experiment. 

In our experiment, we have used the Brazilian test set 
up as indicated in Fig. 1. How the force field varies from 
the point of loading is entirely dependent on the shape and 



Fig. 10. Snapshot of Brazilian test on an-isotropic shale sam¬ 
ple at SINTEF lab [8] 


the material of the loading platen. For simplicity, we as¬ 
sume that the force decreases linearly from the tip point. 
It may also be possible to design loading plates where the 
gradient of the applied force follows power law or expo¬ 
nential. In our simulation we can easily include this type 
of force gradient along the loading plate. We believe that, 
it will have little effect on the nature of the final frac¬ 
ture pattern - because the main fracture initiates at and 
around the most stressed elements and after initial failure 
the newly created weak zones guide the fracture develop¬ 
ment. Therefore, as long as we have a force distribution 
on the system having a stress-tip and symmetrical gradi¬ 
ent (decreasing) at both sides from the tip, we can expect 
similar fracture pattern. 

For breaking thresholds distributed uniformly and with 
Weibull distribution, we observe that the critical strength 
cr c (L ) of a system of size L approaches the asymptotic val¬ 
ues 0.14 and 0.24 respectively. For all the results obtained 
from our model, we have taken the strength decrement 
factor a = 0.7 and the cut-off length of the applied force 
L C = L/ 12. 

The stress-strain curve is very sensitive to the value 
of a. In Fig. [5] the stress-strain behavior of a system with 
L = 64 has been shown for three different values of a. 
Lowering the value of a the more and more narrow, ver¬ 
tically directed and localized the damage zone become, 
shown in Fig. [6l As an effect the strength of the system 
will decrease. For a = 1 the behavior of this model is 
same as the mean field version of the Fiber Bundle model 
mm, the elements break one by one with increasing 
sequence of their breaking thresholds. 

Our statistical model is very simple in the sense -there 
are only two tuning parameters: the strength decrement 
factor a and the slope of the applied force/ load. Ac¬ 
cording to the classical theory -stress fields are enhanced 
around the failure/ defect zones -which in turn weaken 
the elements in those areas. This fact has been captured 
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in our model through the factor a. The second parameter 
can be treated as an initial boundary condition which de¬ 
pends on the geometrical arrangements for the Brazilian 
test. 

Recently Simpson et. al. has done a very interesting 
Brazilian test [SJ on anisotropic rock sample (shale) and 
has monitored fracture development through high speed 
camera (See Fig. fT0|) . Can we apply our model to sim¬ 
ulate/explain such Brazilian test on highly anisotropic 
rocks na? The answer is yes. We can easily create weak 
bedding planes inside the model rock mass and we ex¬ 
pect that in our model when the growing-fracture meets a 
weak plane - it must follow the direction of the weak plane 
unless it gets arrested by the presence of some strong el¬ 
ements. However we have to create the 3D version of our 
model first and then it needs a proper calibration against 
lab-experiments -which is our ongoing research activity. 

The research work in this paper is a part of the sci¬ 
entific collaboration in the INDNOR project 217413/E20 
funded by the Research Council of Norway. S. K. is thank¬ 
ful for financial support through the project grant to visit 
SINTEF Petroleum Research. We thank S. S. Manna for 
useful discussions. S.P. acknowledges partial support from 
SINTEF Petroleum Research through internal SIP project 
(7020606). 
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